rm(list=ls()) 

library(arm)
library(VGAM)

## get estimates from Bayesx output

a.results <- read.table("spatial_f_district_spatial.res",header=T)

b.names <- c("73rd","74th","75th","76th","77th","78th","79th","80th")

b.names.alt <- t(matrix(b.names,3,8))

b.names.alt <- c(matrix(b.names.alt,24,1),"Democrat","Labor Comm","Urban")

b.fe.results <- read.table("spatial_FixedEffects1.res",header=T)

dta <- read.table("../../rollcall.logit.inter.all.txt",header=T)

a <- read.table("spatial_f_district_spatial_sample.raw",header=T)[,-1]

b.union <- read.table("spatial_unionpop_f_district_spatial_sample.raw",header=T)[,-1]

b.aapct <- read.table("spatial_aapct_f_district_spatial_sample.raw",header=T)[,-1]

b.urbanpct <- read.table("spatial_urbanpct_f_district_spatial_sample.raw",header=T)[,-1]

b.fe <- read.table("spatial_FixedEffects1_sample.raw",header=T)

congnum <- seq(73,80,1)

lb <- .05
ub <- .95

regression.plots <- function (ests,varname,regionname) {

lbub <- apply(ests,2,quantile,(probs=c(lb,ub)))

b.range <- range(lbub[1,],lbub[2,])+range(lbub[1,],lbub[2,])*.01

plot(congnum,apply(ests,2,median),ylim=b.range,pch=20,cex.axis=.8,xlab="Congress",ylab="Estimates",main=list(varname,cex=.9),sub=regionname)

segments(congnum,lbub[1,],congnum,lbub[2,])

lines (c(72,81), c(0,0),lty=3)

cat(regionname,varname,"\n")
print(cbind(congnum,apply(ests,2,median),lbub[1,],lbub[2,]))

}

pdf (paste(file="votes_regression_figs.pdf"))

par(mfrow=c(4,3))

regression.plots(a[,1:8],"Region-Period Effect","Deep South")
regression.plots(a[,9:16],"Region-Period Effect","Border South")
regression.plots(a[,17:24],"Region-Period Effect","Non-South")

regression.plots(b.union[,1:8],"Union Pct","Deep South")
regression.plots(b.union[,9:16],"Union Pct","Border South")
regression.plots(b.union[,17:24],"Union Pct","Non-South")

regression.plots(b.aapct[,1:8],"African-American Pct","Deep South")
regression.plots(b.aapct[,9:16],"African-American Pct","Border South")
regression.plots(b.aapct[,17:24],"African-American Pct","Non-South")

regression.plots(b.urbanpct[,1:8],"Urban Pct","Deep South")
regression.plots(b.urbanpct[,9:16],"Urban Pct","Border South")
regression.plots(b.urbanpct[,17:24],"Urban Pct","Non-South")

graphics.off()

## compute marginal effects

## b_1 constant 
## b_2 dem      
## b_3 laborcomm

int <- b.fe$b_1
b.dem <- b.fe$b_2
b.laborcomm <- b.fe$b_3

n.sims <- nrow(a)

n.period <- ncol(a)/3

p.ds1 <- array(NA,c(n.sims,n.period))
p.ds2 <- array(NA,c(n.sims,n.period))
p.aa.ds2 <- array(NA,c(n.sims,n.period))
p.un.aa.ds2 <- array(NA,c(n.sims,n.period))
p.aa.urb.ds2 <- array(NA,c(n.sims,n.period))
p.bs1 <- array(NA,c(n.sims,n.period))
p.bs2 <- array(NA,c(n.sims,n.period))
p.aa.bs2 <- array(NA,c(n.sims,n.period))
p.un.aa.bs2 <- array(NA,c(n.sims,n.period))
p.aa.urb.bs2 <- array(NA,c(n.sims,n.period))
p.ns.dem1 <- array(NA,c(n.sims,n.period))
p.ns.dem2 <- array(NA,c(n.sims,n.period))
p.un.aa.ns.dem2 <- array(NA,c(n.sims,n.period))
p.aa.ns.dem2 <- array(NA,c(n.sims,n.period))
p.aa.urb.ns.dem2 <- array(NA,c(n.sims,n.period))
p.ns.rep1 <- array(NA,c(n.sims,n.period))
p.ns.rep2 <- array(NA,c(n.sims,n.period))

for (j in 1:n.period){

  aapct.med <- median(dta$aapct[dta$period==j & dta$region==1]) 
  urbanpct.med <- median(dta$urbanpct[dta$period==j & dta$region==1])
  unionpop.med <- median(dta$unionpop[dta$period==j & dta$region==1])
  
  aapct.sd <- sd(dta$aapct[dta$period==j & dta$region==1]) 
  urbanpct.sd <- sd(dta$urbanpct[dta$period==j & dta$region==1])
  unionpop.sd <- sd(dta$unionpop[dta$period==j & dta$region==1])
  
    p.ds1[,j] <- probit(int + a[,j] + b.dem +
    b.urbanpct[,j]*urbanpct.med +
    b.union[,j]*unionpop.med +
    b.aapct[,j]*aapct.med, inverse=T)

    p.ds2[,j] <- probit(int + a[,j] + b.dem +
    b.aapct[,j]*aapct.med +
    b.urbanpct[,j]*urbanpct.med +
    b.union[,j]*(unionpop.med + unionpop.sd), inverse=T)

    p.aa.ds2[,j] <- probit(int + a[,j] + b.dem +
    b.aapct[,j]*(aapct.med + aapct.sd) +
    b.urbanpct[,j]*urbanpct.med +
    b.union[,j]*(unionpop.med), inverse=T)

    p.un.aa.ds2[,j] <- probit(int + a[,j] + b.dem +
    b.aapct[,j]*(aapct.med + aapct.sd) +
    b.urbanpct[,j]*urbanpct.med +
    b.union[,j]*(unionpop.med + unionpop.sd), inverse=T)

    p.aa.urb.ds2[,j] <- probit(int + a[,j] + b.dem +
    b.aapct[,j]*(aapct.med + aapct.sd) +
    b.urbanpct[,j]*(urbanpct.med + urbanpct.sd) +
    b.union[,j]*(unionpop.med + unionpop.sd), inverse=T)

}
    
for (j in 1:n.period){

  aapct.med <- median(dta$aapct[dta$period==j & dta$region==2]) 
  urbanpct.med <- median(dta$urbanpct[dta$period==j & dta$region==2])
  unionpop.med <- median(dta$unionpop[dta$period==j & dta$region==2])
  
  aapct.sd <- sd(dta$aapct[dta$period==j & dta$region==2]) 
  urbanpct.sd <- sd(dta$urbanpct[dta$period==j & dta$region==2])
  unionpop.sd <- sd(dta$unionpop[dta$period==j & dta$region==2])

  k <- n.period+8
  
    p.bs1[,j] <- probit(int + a[,k] + b.dem +
    b.aapct[,k]*aapct.med +
    b.urbanpct[,k]*urbanpct.med +
    b.union[,k]*unionpop.med, inverse=T)

    p.bs2[,j] <- probit(int + a[,k] + b.dem +
    b.urbanpct[,k]*urbanpct.med +
    b.aapct[,k]*aapct.med +
    b.union[,k]*(unionpop.med + unionpop.sd), inverse=T)

    p.aa.bs2[,j] <- probit(int + a[,k] + b.dem +
    b.urbanpct[,k]*urbanpct.med +
    b.aapct[,k]*(aapct.med + aapct.sd) +
    b.union[,k]*unionpop.med, inverse=T)

    p.un.aa.bs2[,j] <- probit(int + a[,k] + b.dem +
    b.urbanpct[,k]*urbanpct.med +
    b.aapct[,k]*(aapct.med + aapct.sd) +
    b.union[,k]*(unionpop.med + unionpop.sd), inverse=T)

    p.aa.urb.bs2[,j] <- probit(int + a[,k] + b.dem +
    b.aapct[,k]*(aapct.med + aapct.sd) +
    b.urbanpct[,k]*(urbanpct.med + urbanpct.sd) +
    b.union[,k]*(unionpop.med + unionpop.sd), inverse=T)
  
}

  for (j in 1:n.period){ 

  aapct.med <- median(dta$aapct[dta$period==j & dta$region==3]) 
  urbanpct.med <- median(dta$urbanpct[dta$period==j & dta$region==3])
  unionpop.med <- median(dta$unionpop[dta$period==j & dta$region==3])
  
  aapct.sd <- sd(dta$aapct[dta$period==j & dta$region==3]) 
  urbanpct.sd <- sd(dta$urbanpct[dta$period==j & dta$region==3])
  unionpop.sd <- sd(dta$unionpop[dta$period==j & dta$region==3])

  k <- n.period+16
  
    p.ns.dem1[,j] <- probit(int + a[,k] + b.dem +
    b.aapct[,k]*aapct.med +
    b.urbanpct[,k]*urbanpct.med +
    b.union[,k]*unionpop.med, inverse=T)

    p.ns.dem2[,j] <- probit(int + a[,k] + b.dem +
    b.aapct[,k]*aapct.med +
    b.urbanpct[,k]*urbanpct.med +
    b.union[,k]*(unionpop.med + unionpop.sd), inverse=T)

    p.aa.ns.dem2[,j] <- probit(int + a[,k] + b.dem +
    b.aapct[,k]*(aapct.med + aapct.sd) +
    b.urbanpct[,k]*urbanpct.med +
    b.union[,k]*unionpop.med, inverse=T)

    p.un.aa.ns.dem2[,j] <- probit(int + a[,k] + b.dem +
    b.aapct[,k]*(aapct.med + aapct.sd) +
    b.urbanpct[,k]*urbanpct.med +
    b.union[,k]*(unionpop.med + unionpop.sd), inverse=T)

    p.aa.urb.ns.dem2[,j] <- probit(int + a[,k] + b.dem +
    b.aapct[,k]*(aapct.med + aapct.sd) +
    b.urbanpct[,k]*(urbanpct.med + urbanpct.sd) +
    b.union[,k]*(unionpop.med + unionpop.sd), inverse=T)

    p.ns.rep1[,j] <- probit(int + a[,k] + 
    b.aapct[,k]*aapct.med +
    b.urbanpct[,k]*urbanpct.med +
    b.union[,k]*unionpop.med, inverse=T)

    p.ns.rep2[,j] <- probit(int + a[,k] + 
    b.aapct[,k]*aapct.med +
    b.urbanpct[,k]*urbanpct.med +
    b.union[,k]*(unionpop.med + unionpop.sd), inverse=T)

  
}


p.diff.ds <- p.ds2 - p.ds1
p.diff.bs <- p.bs2 - p.bs1
p.diff.dem.ns <- p.ns.dem2 - p.ns.dem1
p.diff.rep.ns <- p.ns.rep2 - p.ns.rep1

p.diff.aa.ds <- p.aa.ds2 - p.ds1
p.diff.aa.bs <- p.aa.bs2 - p.bs1
p.diff.aa.dem.ns <- p.aa.ns.dem2 - p.ns.dem1

p.diff.un.aa.ds <- p.un.aa.ds2 - p.ds1
p.diff.un.aa.bs <- p.un.aa.bs2 - p.bs1
p.diff.un.aa.dem.ns <- p.un.aa.ns.dem2 - p.ns.dem1

p.diff.aa.urb.ds <- p.aa.urb.ds2 - p.ds1
p.diff.aa.urb.bs <- p.aa.urb.bs2 - p.bs1
p.diff.aa.urb.dem.ns <- p.aa.urb.ns.dem2 - p.ns.dem1

p.dd.bs <- p.diff.aa.urb.bs - p.diff.bs
p.dd.ds <- p.diff.aa.urb.ds - p.diff.ds
p.dd.dem.ns <- p.diff.aa.urb.dem.ns - p.diff.dem.ns

me.plots <- function (ests,varname,regionname) {

lbub <- apply(ests,2,quantile,(probs=c(lb,ub)))

plot(congnum,apply(ests,2,median),ylim=b.range,pch=20,cex.axis=.8,xlab="Congress",ylab=list("Difference in Predicted Probabilities",cex=.9),main=list(varname,cex=1),sub=regionname)

segments(congnum,lbub[1,],congnum,lbub[2,])

lines (c(72,81), c(0,0),lty=3)

cat(regionname,varname,"\n")
print(cbind(congnum,apply(ests,2,median),lbub[1,],lbub[2,]))


}

pdf (paste(file="me_figs.pdf"))

par(mfrow=c(2,2))

b.range <- range(cbind(p.diff.ds,p.diff.bs,p.diff.dem.ns,p.diff.rep.ns))

me.plots(p.diff.ds,"Increase in Union Membership","Deep South")
me.plots(p.diff.bs,"Increase in Union Membership","Border South")
me.plots(p.diff.dem.ns,"Increase in Union Membership","Non-South, Democrats")
me.plots(p.diff.rep.ns,"Increase in Union Membership","Non-South, Republicans")

pdf (paste(file="me_aa_figs.pdf"))

par(mfrow=c(4,3))

b.range <- range(cbind(p.diff.ds,p.diff.bs,p.diff.dem.ns))
                 
me.plots(p.diff.ds,"Increase in Union Membership","Deep South")
me.plots(p.diff.bs,"Increase in Union Membership","Border South")
me.plots(p.diff.dem.ns,"Increase in Union Membership","Non-South, Democrats")

b.range <- range(apply(cbind(p.diff.aa.ds,p.diff.aa.bs,p.diff.aa.dem.ns),2,quantile,(probs=c(lb,ub))))
                 
me.plots(p.diff.aa.ds,"Increase in African-American Pct","Deep South")
me.plots(p.diff.aa.bs,"Increase in African-American Pct","Border South")
me.plots(p.diff.aa.dem.ns,"Increase in African-American Pct","Non-South, Democrats")

b.range <- range(apply(cbind(p.diff.un.aa.ds,p.diff.un.aa.bs,p.diff.un.aa.dem.ns),2,quantile,(probs=c(lb,ub))))
                 
me.plots(p.diff.un.aa.ds,"Increase in Union Membership \n & African-American Pct","Deep South")
me.plots(p.diff.un.aa.bs,"Increase in Union Membership \n & African-American Pct","Border South")
me.plots(p.diff.un.aa.dem.ns,"Increase in Union Membership \n & African-American Pct","Non-South, Democrats")

b.range <- range(apply(cbind(p.diff.aa.urb.ds,p.diff.aa.urb.bs,p.diff.aa.urb.dem.ns),2,quantile,(probs=c(lb,ub))))

me.plots(p.diff.aa.urb.ds,"Increase in Union Membership, \n African-American & Urban Pct","Deep South")
me.plots(p.diff.aa.urb.bs,"Increase in Union Membership, \n African-American & Urban Pct","Border South")
me.plots(p.diff.aa.urb.dem.ns,"Increase in Union Membership, \n African-American & Urban Pct","Non-South, Democrats")

pdf (paste(file="me_dd_figs.pdf"))

par(mfrow=c(1,3))

b.range <- range(cbind(p.dd.ds,p.dd.bs,p.dd.dem.ns))
                 
    me.plots(p.dd.ds,"Difference in difference","Deep South")
    me.plots(p.dd.bs,"Difference in difference","Border South")
    me.plots(p.dd.dem.ns,"Difference in difference","Non-South, Democrats")


graphics.off()

a.names <- c("73rd,Deep South","73rd,Border South","73rd,Non-South","74th,Deep South","74th,Border South","74th,Non-South","75th,Deep South","75th,Border South","75th,Non-South","76th,Deep South","76th,Border South","76th,Non-South","77th,Deep South","77th,Border South","77th,Non-South","78th,Deep South","78th,Border South","78th,Non-South","79th,Deep South","79th,Border South","79th,Non-South","80th,Deep South","80th,Border South","80th,Non-South")

a.names.alt <- t(matrix(a.names,3,8))

par(mfrow=c(1,2))

rnames <- matrix(a.names.alt,1,24)

texfile <- "rc_probits_mrf.tex"

roundval <- 3

## make table with results

  cat("\\begin{table}\\centering \\small \\ssp \n ",file=texfile,append=FALSE)
  cat("\\begin{threeparttable}\n ",file=texfile,append=TRUE)
  cat("\\caption{Results from analysis of labor roll call votes in the Senate using Markov field random priors (73rd--80th Congresses)} \n \\label{tab:rc_probits_mrf} \n",file=texfile,append=TRUE)
  cat("\\begin{tabular}{ld{-1}d{-1}d{-1}d{-1}} \n \\toprule \n \\midrule \n",file=texfile,append=TRUE)
  cat("& \\multicolumn{1}{c}{Region-Period Effect} & \\multicolumn{1}{c}{Union} & \\multicolumn{1}{c}{AA\\%} & \\multicolumn{1}{c}{Urban \\%} \\\\ \n",file=texfile,append=TRUE)

  names.ds <- rnames[,1]

## compute median and credible interval and see if latter include 0

## loop through each region and each coefficient

ests <- cbind(apply(a,2,median),apply(b.union,2,median),apply(b.aapct,2,median),apply(b.urbanpct,2,median))

a.sig <- ifelse(apply(a,2,quantile,(probs=c(lb,ub)))[1,] <= 0 & apply(a,2,quantile,(probs=c(lb,ub)))[2,] >= 0,"","^{*}")

b.un.sig <- ifelse(apply(b.union,2,quantile,(probs=c(lb,ub)))[1,] <= 0 & apply(b.union,2,quantile,(probs=c(lb,ub)))[2,] >= 0,"","^{*}")

b.aa.sig <- ifelse(apply(b.aapct,2,quantile,(probs=c(lb,ub)))[1,] <= 0 & apply(b.aapct,2,quantile,(probs=c(lb,ub)))[2,] >= 0,"","^{*}")

b.urb.sig <- ifelse(apply(b.urbanpct,2,quantile,(probs=c(lb,ub)))[1,] <= 0 & apply(b.urbanpct,2,quantile,(probs=c(lb,ub)))[2,] >= 0,"","^{*}")

sigmat <- cbind(a.sig,b.un.sig,b.aa.sig,b.urb.sig)

  for (t in 1:length(rnames)){
    if (t == 9 | t == 17) cat("\n \\midrule \n",file=texfile,append=TRUE)
    cat(rnames[t],file=texfile,append=TRUE)
    for (k in 1:ncol(ests)){
      cat (paste(" & ", round(ests[t,k],roundval),sigmat[t,k]),file=texfile,append=TRUE)
    }
    cat (" \\\\\n ",file=texfile,append=TRUE)
  }

cat("\n \\midrule \\midrule \n",file=texfile,append=TRUE)

cat(paste("Democrat &", round(median(b.dem),roundval),ifelse(quantile(b.dem,(probs=c(lb,ub)))[1] <= 0 & apply(a,2,quantile,(probs=c(lb,ub)))[2] >= 0,"","^{*}")," \\\\ \n"),file=texfile,append=TRUE)

cat(paste("Labor Committe &", round(median(b.laborcomm),roundval),ifelse(quantile(b.dem,(probs=c(lb,ub)))[1] <= 0 & apply(a,2,quantile,(probs=c(lb,ub)))[2] >= 0,"","^{*}")," \\\\ \n"),file=texfile,append=TRUE)

cat("\\bottomrule \n \\end{tabular}\n \\begin{tablenotes} \n \\item \n {\\em Notes}: \\ssp $^{*} p \\leq .1$. Estimation using MCMC implemented in {\tt BayesX}, 100,000 iterations (first 10,000 discarded). \\textit{Union} indicates union density, \\textit{AA} indicates percent African-American, \\textit{Urban} indicates percent urbanization, and \\textit{Labor Committee} and \\textit{Democrat} refer to dummy variables indicating membership on the labor committee and in Democratic party, respectively. \n \\end{tablenotes} \n \\end{threeparttable} \n \\hspace{\\fill} \n \\end{table}\n \n \n",file=texfile,append=TRUE)

  
